rm(list=ls())
library(maptools)
library(sp)
library(plyr)
library(rgdal)
library(lme4)
library(ggplot2)
library(doBy)
library(date)
library(foreign)
library(Hmisc)
library(memisc)
library(mgcv)
library(rms)
library(lmtest)
library(sandwich)
library(rdrobust)
library(apsrtable)


wd<-"YOUR WORKING DIRECTORY HERE"
setwd(wd)

##define function to estimate local effect at various bandwidths
source(paste(wd,"functions/","local_discont_function.R",sep=""))

##define function to generate regression output on large models
source("functions/full.models.R")

##define function to plot loess pre/post cutoff date
source("functions/loess.time.plot.R")

##define function to retrieve larger SE between conventional and HAC
source("functions/larger.se.R")

##define 5-fold cross validation function
source(paste(wd,"functions/","cross_val5.R",sep=""))

##define function for clustered SEs
source(paste(wd,"functions/","vcovCluster.R",sep=""))


#### LOAD DATA SETS
load(file="data/sf_2.Rdata") ##all stops 2008-2015, "s"
s<-s[order(s$dayno),]
dim(s)
table(is.na(s$lag.hit))
range(s$dayno)
range(as.numeric(s$dayno))
identical(s$dayno, s$dayno[order(s$dayno)])
range(as.numeric(s$dayno))

load("data/sw_2.Rdata")##weapon stops 2008-2015, "sw"
dim(sw)
sw<-sw[order(sw$dayno),]
table(is.na(sw$lag.hit))
range(sw$dayno)
range(as.numeric(sw$dayno))
identical(sw$dayno, sw$dayno[order(sw$dayno)])
#write.dta(sw, file="data/sw.dta")

load(file="data/sb_2.Rdata")##all stops 2008-2015 with block group data, "sb"
dim(sb)
sb<-sb[order(sb$dayno),]
table(is.na(sb$lag.hit))
range(sb$dayno)
range(as.numeric(sb$dayno))
identical(sb$dayno, sb$dayno[order(sb$dayno)])


dates<-c("1Jan2008","1Jan2010","1Jan2012","1Jan2014","1Jan2016")


##define date of intervention
memo<-mdy.date(month=3, year=2013, day=5)#define intervention day, March 5, 2013


##Note the data are ordered by dayno, which is required for the use of vcovHAC

table(sw$distance[sw$dayno==(memo-2)])##days before the intervention are negative
table(sw$distance[sw$dayno==memo])##day of intervention is zero
table(sw$distance[sw$dayno==(memo+2)])##days after intervention are positive



##############
############# ANALYSIS
##############
###############


###############
###############
#### Estimates quoted in text that don't appear in a table or figure:
###############
###############


##Share of stops in pre-treatment period that resulted in arrest or summons
table(s$arrested)
table(s$summons.issued)
table(I(s$arrested==T|s$summons.issued==T))
mean(I(s$arrested==T|s$summons.issued==T), na.rm=T)
mean(I(s$arrested[s$year<2013]==T|s$summons.issued[s$year<2013]==T), na.rm=T)


##Share of stops that are of nonwhite suspects
table(s$race2)
mean(I(s$race2=="white"), na.rm=T)
mean(I(s$race2[s$year<2013]=="white"), na.rm=T)


##How many weapon stops in the data? (Discussed in "Results" section)
sum(table(sw$weapon.hit))##826,573
##What % of all stops are weapon stops?
mean(I(s$suspected.crime=="cpw"), na.rm=T)##25.9%


##how many weapons found by year?
r<-summaryBy(found.weapon~year, data=s, FUN=sum, na.rm=T)
r
##how many guns by year?
r<-summaryBy(found.gun~year, data=s, FUN=sum, na.rm=T)
r
##how many knives/other
r<-summaryBy(found.knife+found.other~year, data=s, FUN=sum, na.rm=T)
r
##means
r<-summaryBy(found.gun +found.knife+found.other~year, data=s, FUN=mean, na.rm=T)
r
r<-summaryBy(found.gun +found.knife+found.other~year, data=sw[sw$weapon.hit==1 & sw$dayno<memo,], FUN=mean, na.rm=T)
r

mean(sw[sw$weapon.hit==1 & sw$dayno<memo,"found.gun"],na.rm=T)

############################
############################ What is the mean hit rate during the pre-treatment period?
############################
############################

mean(sw$weapon.hit[sw$dayno<memo])




############################
############################
############################
############################ TABLE 1
############################
############################
############################



##TABLE 1
## Estimate discontinuity at moment of intervention using all weapon stops, 2008-2015

## Executing the following line of code will estimate the var-covar matrix using HAC standard errors for all models and save it in the "output" folder. This is very time consuming, so the estimated matrix is already stored in the "output" folder.
#source("batch/weapon_models_batch.R")

##To estimate standard (homoscedastic) SEs and save on estimation time, set "std" equal to T in the full.data.models function
res<-full.data.models(data=sw, dv="weapon.hit", short.data=F,  std=T)
res

##To recover the larger of HAC and conventional SEs, supply the path to the file containing the robust SEs in the "vc" parameter
res<-full.data.models(data=sw, dv="weapon.hit", short.data=F,  vc="output/weapon_HAC_vc.Rdata", std=F)
##print LaTex table
res[["table"]]##this matches table 1 in the paper




############################
############################
############################
############################ FIGURE 1
############################
############################
############################

#DAILY WEAPON HIT RATE OVER TIME

##compute daily hit rate in the data pre/post intervention using loess estimator
out.hit<-loess.time.series(data=sw, dv="weapon.hit",group="dayno", fun=mean, cut=memo)
##compute daily numerator and denominator of hit rate over time
sw$count<-1
out.denom<-loess.time.series(data=sw, dv="count",group="dayno", fun=sum, cut=memo)
out.num<-loess.time.series(data=sw, dv="weapon.hit",group="dayno", fun=sum, cut=memo)

#generate figure

file.name<-paste(wd, "output/hit.rate.time_loess_panel_a.pdf", sep="")
pdf(file.name)
#par(mar=c(4,4,4,4))

means.all.day<-out.hit$means.all.day
dd1<-out.hit$dd1
dd2<-out.hit$dd2
fit1<-out.hit$lower$fit1
fit2<-out.hit$upper$fit2
i1<-out.hit$lower$i1
i2<-out.hit$upper$i2


upper<-quantile(out.hit$means.all.day$weapon.hit.fun, probs=c(.99))

plot(means.all.day$dayno, means.all.day[,2], pch=19, col="grey", cex=.3, xlim=c(min(means.all.day$dayno, na.rm=T), max(means.all.day$dayno, na.rm=T)), xlab="", ylab="", main="Hit Rate Increases Markedly\nWhen Procedural Reform Introduced", axes=F, cex.main=1.3, ylim=c(0,upper))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=T, las=2, col.axis="black"
)
mtext(side=2, "Daily Weapon Recovery Rate", line=3)
mtext(side=1, "Date", line=3, cex=.8)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
text( x=mdy.date(month=3, year=2013, day=5)-440, y=.24,labels="Day of\nProcedural\nChange" )
lines(dd1$dayno[i1], fit1[i1], col="black", lwd=2)
lines(dd2$dayno[i2], fit2[i2], col="black", lwd=2)
arrows(mdy.date(month=3, year=2013, day=5)-400, .22,mdy.date(month=3, year=2013, day=5), .22, length=.13 )

dev.off()

file.name<-paste(wd, "output/hit.rate.time_loess_panel_b.pdf", sep="")
pdf(file.name, width=8, height=4)

layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(1,1), heights=c(1,1))

##global numerator
means.all.day<-out.num$means.all.day
dd1<-out.num$dd1
dd2<-out.num$dd2
fit1<-out.num$lower$fit1
fit2<-out.num$upper$fit2
i1<-out.num$lower$i1
i2<-out.num$upper$i2
plot(means.all.day$dayno, means.all.day[,2], pch=19, col="grey", cex=.3,  axes=F, xlab="Date",ylab="# Daily Stops Producing a Weapon", ylim=c(0,max(means.all.day[,2])), lty=1,lwd=2, main="Stops Producing Weapons \n(Numerator)", xlim=c(min(means.all.day$dayno), max(means.all.day$dayno) ) )
lines(dd1$dayno[i1], fit1[i1], col="black", lwd=2, lty=1)
lines(dd2$dayno[i2], fit2[i2], col="black", lwd=2, lty=1)
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, at=seq(0, 50, by=10), las=2, cex.axis=.8)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo



##global denominator
means.all.day<-out.denom$means.all.day
dd1<-out.denom$dd1
dd2<-out.denom$dd2
fit1<-out.denom$lower$fit1
fit2<-out.denom$upper$fit2
i1<-out.denom$lower$i1
i2<-out.denom$upper$i2
plot(means.all.day$dayno, means.all.day[,2], pch=19, col="grey", cex=.3,  axes=F, xlab="Date",ylab="# Daily Weapon Stops", ylim=c(0,max(means.all.day[,2])), lty=1,lwd=2, main="Weapon Stops\n(Denominator)", xlim=c(min(means.all.day$dayno), max(means.all.day$dayno) ) )
lines(dd1$dayno[i1], fit1[i1], col="black", lwd=2, lty=1)
lines(dd2$dayno[i2], fit2[i2], col="black", lwd=2, lty=1)
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, at=seq(0, 1000, by=100), las=2, cex.axis=.8)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo


dev.off()




############################
############################
############################
############################ FIGURE 2
############################
############################
############################


##LOCAL ESTIMATES OF DISCONTINUITY
stops.new.short<-sw[sw$dayno>=(memo-30) & sw$dayno<(memo+30), ]

identical(stops.new.short$dayno,stops.new.short$dayno[order(stops.new.short$dayno)] )

dim(stops.new.short)
length(unique(stops.new.short$dayno)) ##30 days on either side of treatment

# out<-models(data=stops.new.short, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols",int="no",cluster=F)###


# save(out, file="output/figure3_data.Rdata")

load("output/figure3_data.Rdata")

##view results for one model
out[["m1"]]

m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-1:30

##All models
file.name<-"output/bandwidth.pdf"
pdf( file.name, width=6, height=6)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1))
plot( x.axis, m5$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m1$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m2$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m3$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m4$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic 
+ Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m6$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)


dev.off()



############################
############################
############################TABLE 2
############################
############################


sw2<-na.omit(sw[,c("period","distance","lag.hit","dayno","weapon.hit","weekday","year","month")])
sw2$count<-1
head(sw2)

##total weapon stops
stops<-summaryBy(count~dayno, data=sw2, FUN=sum, na.rm=T)
##total weapon stops producing a weapon
hits<-summaryBy(weapon.hit~dayno, data=sw2, FUN=sum, na.rm=T)

sw2$period2<-NA
sw2$period2[sw2$period=="pre-memo"]<-0
sw2$period2[sw2$period=="post-memo"]<-1
table(sw2$period2)

temp<-sw2[,c("period2","weekday","year","month","lag.hit","distance","dayno")]
temp<-unique(temp)
dim(temp)

sw3<-merge(hits, stops, by="dayno",all.x=T)
sw3<-merge(sw3, temp, by="dayno",all.x=T)
colnames(sw3)<-c("dayno","weapon.hit.sum","count.sum","period","weekday","year","month","lag.hit","distance")
head(sw3)
dim(sw3)

sw3$period[sw3$period==0]<-"pre-memo"
sw3$period[sw3$period==1]<-"post-memo"
sw3$period<-factor(sw3$period, levels=c("pre-memo","post-memo"))
table(sw3$period)

head(sw3)

sw3<-sw3[order(sw3$dayno),]
res.num<-full.data.models(data=sw3, dv="weapon.hit.sum", short.data=T, short.data.cut=100,  vc=NULL, gen.vc=T)
##print LaTex table
res.num[["table"]]

res.denom<-full.data.models(data=sw3, dv="count.sum", short.data=T, short.data.cut=100,  vc=NULL, gen.vc=T)
##print LaTex table
res.denom[["table"]]








#######################
#######################
#######################
####################### FIGURE 3 and 4 , HOMICIDES
#######################
#######################
#######################
###nyc annual murders

##source: FBI UCR
d<-read.csv("data/nyc_violent_crime.csv",header=T)
head(d)

year<-as.numeric(as.character(d$Year))
year[length(year)+1]<-2015
year

##manually add data for extra years
mur<-as.numeric(as.character(d$Murder.and.nonnegligent.Manslaughter))
mur[length(mur)+1]<-352 #2015 source: http://www.nyc.gov/html/nypd/downloads/pdf/analysis_and_planning/seven_major_felony_offenses_2000_2015.pdf

mur

file.name<-"output/murders_ann.pdf"
pdf(file=file.name, height=6, width=8)

plot(year, mur, type="b", pch=19, main="Annual Murders in NYC", yaxt="n", xlab="Year",ylab="Number of Murders")
axis(2, las=2, at=seq(0, 3000, by=250))
points(year, mur, pch=19)
abline(v=2013, lty=2, col="black")
text( x=2010, y=1750,labels="Year of\nProcedural\nChange" )
arrows(2010, 1600,2013, 1600, length=.13 )

dev.off()



##homicides pre/post

d<-read.csv("data/nypd_weekly_homicides_2011_2013.csv")
head(d)
d$num<-1:nrow(d)
d
treat.location<-d$num[d$date=="3/4/2013-3/10/2013"]
d$treat<-0
d$treat[d$num>= treat.location]<-1
table(d$treat)
d$weeks.since<-(1-treat.location):((nrow(d)-treat.location))
##check coding of distance variable
d$weeks.since[d$num==treat.location-1]#is this -1? yes
d$weeks.since[d$num==treat.location]#is this zero? yes
d$weeks.since[d$num==treat.location+1]#is this 1? yes
##coding correct

d$date<-as.character(d$date)
temp<-strsplit(d$date, "-")
d$start.day<-NA
d$start.month<-NA
d$start.year<-NA

d$stop.day<-NA
d$stop.month<-NA
d$stop.year<-NA

for(i in 1:length(temp)){
	
d$start.month[i]<-strsplit(temp[[i]][1], "/")[[1]][1]
d$start.day[i]<-strsplit(temp[[i]][1], "/")[[1]][2]
d$start.year[i]<-strsplit(temp[[i]][1], "/")[[1]][3]

d$stop.month[i]<-strsplit(temp[[i]][2], "/")[[1]][1]
d$stop.day[i]<-strsplit(temp[[i]][2], "/")[[1]][2]
d$stop.year[i]<-strsplit(temp[[i]][2], "/")[[1]][3]
	
}
d$start.year[d$start.year=="3013"]<-"2013"##fix typo

head(d)



##fit pre-treat loess
model1<-predict(loess(homicides ~ weeks.since, data=d[d$treat==0,], span=.5, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1h<-model1$fit
lb1h<-model1$fit-1.96*model1$s
ub1h<-model1$fit+1.96*model1$s
i1h <- order(d[d$treat==0,"weeks.since"])

model2<-predict(loess(homicides ~ weeks.since, data=d[d$treat==1,], span=.8, model=TRUE, x=TRUE, y=TRUE), se=T)
fit2h<-model2$fit
lb2h<-model2$fit-1.96*model2$s
ub2h<-model2$fit+1.96*model2$s
i2h <- order(d[d$treat==1,"weeks.since"])





#residualize with respect to month indicators
d$my<-(m<-lm(homicides~factor(start.month), data=d))$residuals


##fit pre-treat loess
model1<-predict(loess(my ~ weeks.since, data=d[d$treat==0,], span=.5, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1hh<-model1$fit
lb1hh<-model1$fit-1.96*model1$s
ub1hh<-model1$fit+1.96*model1$s
i1hh <- order(d[d$treat==0,"weeks.since"])


##fit post-treat loess. data is sparser so am using a larger span.
model2<-predict(loess(my ~ weeks.since, data=d[d$treat==1,], span=.8, model=TRUE, x=TRUE, y=TRUE), se=T)
fit2hh<-model2$fit
lb2hh<-model2$fit-1.96*model2$s
ub2hh<-model2$fit+1.96*model2$s
i2hh <- order(d[d$treat==1,"weeks.since"])





file.name<-"output/homicide.pdf"
pdf(file.name, width=12, height=5)

par(mar=c(4,6,2,4), mfrow=c(1,2))
plot(d$weeks.since, d$homicides, col="grey",cex=.6, axes=F,pch=19, ylab="Weekly Homicides",xlab="Weeks From Intervention", main="Weekly Homicides\n(Raw Data)", cex.main=1, cex.axis=.5)
lines(d$weeks.since[d$treat==0][i1h], fit1h[i1h], col="black",lwd=2)
lines(d$weeks.since[d$treat==1][i2h], fit2h[i2h], col="black",lwd=2)
abline(v=0, lty=2, col="black", lwd=1)
axis(1, at=seq(-200, 100, by=50) , cex.axis=1)
axis(2, at=seq(0,500,5),las=2, cex.axis=1)

plot(d$weeks.since, d$my, col="grey",cex=.6, axes=F,pch=19, ylab="Weekly Homicides\n(Controlling for Month)",xlab="Weeks From Intervention", main="Weekly Homicides\n(Controlling for Month)", cex.main=1, cex.axis=.5)
lines(d$weeks.since[d$treat==0][i1h], fit1hh[i1hh], col="black",lwd=2)
lines(d$weeks.since[d$treat==1][i2h], fit2hh[i2hh], col="black",lwd=2)
abline(v=0, lty=2, col="black", lwd=1)
axis(1, at=seq(-200, 100, by=50) , cex.axis=1)
axis(2, at=seq(-500,500,5),las=2, cex.axis=1)


dev.off()





##############################
##############################
##############################
##############################
############################## FIGURE 4: Checks for Reporting Bias
##############################
##############################
##############################
##############################


### UNIFORM ANALYSIS


stops.new.short<-sw[sw$dayno>=(memo-30) & sw$dayno<(memo+30) & sw$officer.uniform==T, ]
# out<-models(data=stops.new.short, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols",int="no", cluster=F)
# save(out, file="output/uniform_data.Rdata")

load("output/uniform_data.Rdata")
##reorganize results for plotting
#relabel
m1a<-out[["m1"]]
m2a<-out[["m2"]]
m3a<-out[["m3"]]
m4a<-out[["m4"]]
m5a<-out[["m5"]]
m6a<-out[["m6"]]



##Change among nonweapon stops

s$cpw<-I(s$suspected.crime=="cpw")
stops.new<-s
stops.new.short<-stops.new[stops.new$dayno>=(memo-30) & stops.new$dayno<(memo+30) & stops.new$cpw==F, ]
# out<-models(data=stops.new.short, N=1:30, dv="found.weapon",memo=memo, model="ols",int="no",cluster=F)
# save(out, file="output/figure3_data_noncpw.Rdata")

load("output/figure3_data_noncpw.Rdata")
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

### change in number of crimes
##make indicator for weapon stop
stops.new<-s

#Estimate discontinuity among non-weapon stops
stops.new$cpw<-0
stops.new$cpw[stops.new$suspected.crime=="cpw"]<-1
stops.new$cpw[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$cpw)

stops.new$notcpw<-0
stops.new$notcpw[stops.new$suspected.crime!="cpw"]<-1
stops.new$notcpw[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$notcpw)

stops.new$robbery<-0
stops.new$robbery[stops.new$suspected.crime=="robbery"]<-1
stops.new$robbery[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$robbery)

stops.new$burglary<-0
stops.new$burglary[stops.new$suspected.crime=="burglary"]<-1
stops.new$burglary[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$burglary)




##top 10 crimes
crimes<- c(    "petit larceny" , "criminal possesion of controlled substance" ,"assault" , "criminal possession of marihuana",    "grand larceny" , "criminal trespass"   ,  "grand larceny auto" , "burglary" ,  "robbery" ,"cpw")

stops.new$other.crime<-as.numeric(I(!(stops.new $suspected.crime%in%crimes) & !is.na(stops.new $suspected.crime)))
table(stops.new $other.crime)
mean(stops.new $other.crime)

s$other.crime<-as.numeric(I(!(s $suspected.crime%in%crimes) & !is.na(s $suspected.crime)))
table(s $other.crime)
mean(s $other.crime)

table(stops.new$suspected.crime)[order(table(stops.new$suspected.crime))]


mean(I(s$suspected.crime%in%crimes))
mean(I(s$suspected.crime[s$dayno<memo]%in%crimes))##over 90% of suspected crimes fall in these categories


crimes<-c(crimes, "other.crime")

length(crimes) ##top 10 crimes plus "other" category


#subset to 100 days pre/post
stops.new.short<-stops.new[stops.new$dayno>=(memo-100) & stops.new$dayno<(memo+100),]
length(unique(stops.new.short$dayno))
stops.new.short$count<-1

diffs<-NA
se<-
lb<-NA
ub<-NA
res<-cbind.data.frame(crime=crimes, pre.sum=NA, post.sum=NA, pct.change=NA,diffs, se, lb, ub)
colnames(res)<-c("crime","pre.sum","post.sum","pct.change","mean","se","lb","ub")
head(res)



for(i in 1:length(crimes)){
	
	if(crimes[i]!="other.crime"){
	temp<-stops.new.short[stops.new.short$suspected.crime==crimes[i],]

	sums<-na.omit(summaryBy(count~dayno, data=temp, FUN=sum, na.rm=T))
	sums$period<-I(sums$dayno>=memo)
	sums$distance<-sums$dayno-memo
	
	res$pre.sum[i]<-sum(sums$count.sum[sums$period==F])
	res$post.sum[i]<-sum(sums$count.sum[sums$period==T])
	##estimate linear model
	m<-lm(count.sum~period+distance+period*distance, data=sums)
	v<-coeftest(m, vcov.=vcovHAC(m))
se1<-v[2,"Std. Error"]
se2<-summary(m)$coefficients[2,"Std. Error"]

if(se1>se2){res$se[i]<-se1}
if(se1<se2){res$se[i]<-se2}

res$mean[i]<-summary(m)$coefficients[2,"Estimate"]
res$lb[i]<-res$mean[i]-1.96*res$se[i]
res$ub[i]<-res$mean[i]+1.96*res$se[i]}

if(crimes[i]=="other.crime"){
	
temp<-stops.new.short[stops.new.short[,crimes[i]]==1,]

sums<-na.omit(summaryBy(count~dayno, data=temp, FUN=sum, na.rm=T))
sums$period<-I(sums$dayno>=memo)
sums$distance<-sums$dayno-memo
	
res$pre.sum[i]<-sum(sums$count.sum[sums$period==F])
res$post.sum[i]<-sum(sums$count.sum[sums$period==T])
	
m<-lm(count.sum~period+distance+period*distance, data=sums)
v<-coeftest(m, vcov.=vcovHAC(m))
se1<-v[2,"Std. Error"]
se2<-summary(m)$coefficients[2,"Std. Error"]

if(se1>se2){res$se[i]<-se1}
if(se1<se2){res$se[i]<-se2}

res$mean[i]<-summary(m)$coefficients[2,"Estimate"]
res$lb[i]<-res$mean[i]-1.96*res$se[i]
res$ub[i]<-res$mean[i]+1.96*res$se[i]
	
}

}

res$pct.change<-(res$mean)/(res$pre.sum/100)
res



##make dummies for suspect attributes

vars<-c("suspect.height","suspect.weight","suspect.age","observation.period")


for(i in 1:length(vars)){
	
	stops.new[,paste(vars[i],"2",sep="")]<-as.numeric(as.character(stops.new[,vars[i]]))
	new<-paste(vars[i],"2",sep="")
	stops.new[,new]<-I(stops.new[, new]>=median(stops.new[, new], na.rm=T))	
}


summary(stops.new$suspect.height2)
summary(stops.new$suspect.weight2)
summary(stops.new$suspect.age2)
summary(stops.new$observation.period2)


res<-res[order(res$mean),]
res

crime.label<-c("robbery","weapon","criminal trespass","grand larceny auto","burglary","grand larceny","marijuana","other crime","assault","petit larceny","controlled substance")


### Logit
### propensity scores

##30-day bandwidth

m<-30

stops.new $post<-NA
stops.new $post[stops.new $dayno>=(memo-m) & stops.new $dayno<memo]<-0
stops.new $post[stops.new $dayno>=memo & stops.new $dayno<(memo+m) ]<-1
table(stops.new $post)



d<-na.omit(stops.new[,c("cpw","inside.outside","precinct", "location.housing", "observation.period2", "stopped.bc.object", "stopped.bc.desc", "stopped.bc.casing", "stopped.bc.lookout", "stopped.bc.clothing", "stopped.bc.drugs", "stopped.bc.furtive", "stopped.bc.violent", "stopped.bc.bulge", "stopped.bc.other",  "additional.proximity",  "additional.associating", "additional.direction", "additional.highcrime", "additional.time", "additional.sights", "additional.other", "suspect.sex", "suspect.race", "suspect.age2", "suspect.height2", "suspect.weight2", "suspect.hair","suspect.build","officer.uniform","weekday","hour","dayno","post")])
dim(d)
d$id<-1:(length(d[,1]))
j<-order(d$id[d$post==0])
k<-order(d$id[d$post==1])


psout <- glm(cpw~ inside.outside + factor(precinct) + factor(location.housing) +as.numeric(observation.period2) + stopped.bc.object   +stopped.bc.desc+         stopped.bc.casing+         stopped.bc.lookout+         stopped.bc.clothing+        stopped.bc.drugs+        
stopped.bc.furtive+         stopped.bc.violent+         stopped.bc.bulge+ stopped.bc.other  +additional.proximity          +additional.associating+ additional.direction+ additional.highcrime     
+ additional.time            +additional.sights +          additional.other   +suspect.sex+    suspect.race              
 +suspect.age2+                      +suspect.height2+            
 +suspect.weight2            +suspect.hair              +                suspect.build +officer.uniform+factor(weekday)+factor(hour),family=binomial,data=d[d$post==0, ]) 
summary(psout)
d$pre.prop<-NA
d$pre.prop[d$post==0 ][j]<-psout$fitted.values[j]



##fit a model using the post-treatment observations
prob2<-predict(psout, newdata=d[d$post==1, ], type="response",se.fit=T)

d$post.prop<-NA
d$post.prop[d$post==1 ][k]<-prob2$fit[k]
d$count<-1
sum(d$count[d$post==0])==nrow(d)-table(is.na(d$pre.prop))[2]
sum(d$count[d$post==1])==nrow(d)-table(is.na(d$post.prop))[2]

table(I(is.na(d$pre.prop) & d$post==0))#should be all false
table(I(is.na(d$post.prop) & d$post==1))#should be all false



########PROPENSITY SCORE ANALYSIS
####USE PRE TREAT PARAMETERS TO PREDICT POST-TREATMENT PROB OF BEING LABELED A WEAPON STOP (among non weapon stops)


x.axis<-1:30

y<-nrow(res):1


file.name<-"output/bandwidth_noncpw_prop_uniform.pdf"
pdf( file.name, width=8, height=8)

layout(matrix(c(1,1,2,2,3,3,3,3,4,5,6,7,8,8,8,8,9,10,11,12), ncol=4, byrow = TRUE), widths=c(rep(1,4)), 
heights=c(1,.2,1,.2,1))

par(mar=c(4,7.4,4,4))
plot(res$mean, y, pch=19, cex=.7, axes=F, xlab="Change in No. of Stops", main="Change in Suspected Crimes\non 3/5/2013, 100-Day Bandwidth", ylab="", xlim=c(min(res$lb),0), cex.main=.95)
abline(v=0, lty=2)
segments(x0=res$lb,y0=y, x1=res$ub, y1=y)
axis(1, seq(-200, 1, by=.5)*100, cex.axis=.7)
axis(2, at=y, crime.label , las=2, cex.axis=.8)

plot(density(na.omit(d$pre.prop[d$cpw==0 & d$post==0])), lty=1, ylim=c(0,7.5),main="Prob. of Being Labeled a Weapon Stop\n30-Day Bandwidth", xlab="Propensity Score", yaxt="n")
lines(density(na.omit(d$post.prop[d$cpw==0 &d$post==1])), lty=2, col="red")
legend("topright",legend=c("Pre-Treatment","Post-Treatment"), lty=c(1,2),col=c("black","red"))
abline(v=mean(d$pre.prop[d$cpw==0 & d$post==0], na.rm=T))
abline(v=mean(d$post.prop[d$cpw==0 & d$post==1], na.rm=T), lty=2, col="red")
axis(2,at=seq(0,8,2),las=2)

##insert text
par(mar = c(0,0,1,0)) 
plot(1,1,type = "n",frame.plot = FALSE,axes = FALSE)
u <- par("usr")
text(1,u[4],labels = "Discontinuity Among Non-Weapon Stops",pos = 1, cex=1.5)


par(mar=c(4,4,4,1))
plot( x.axis, m5$coef*100, ylim=c(-3,2), pch=19, xlab="Bandwidth (days)",ylab="",main="Difference in Means", cex=.5, cex.axis=.8, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
axis(2, at=seq(-10,10,by=1), las=2, cex.axis=.7)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)

plot( x.axis, m1$coef*100, ylim=c(-3,2), pch=19, xlab="Bandwidth (days)",ylab="",main="Linear", cex=.5, cex.axis=.8, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
axis(2, at=seq(-10,10,by=1), las=2, cex.axis=.7)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)


plot( x.axis, m2$coef*100, ylim=c(-3,2), pch=19, xlab="Bandwidth (days)",ylab="",main="Quadratic", cex=.5, cex.axis=.7, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
axis(2, at=seq(-10,10,by=1), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)

plot( x.axis, m6$coef*100, ylim=c(-3,2), pch=19, xlab="Bandwidth (days)",ylab="",main="Cubic", cex=.5, cex.axis=.8, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
axis(2, at=seq(-10,10,by=1), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)



par(mar = c(0,0,1,0)) 
plot(1,1,type = "n",frame.plot = FALSE,axes = FALSE)
u <- par("usr")
text(1,u[4],labels = "Discontinuity Among Stops Made By Officers in Uniform",pos = 1, cex=1.5)



par(mar=c(4,4,4,2))
##uniform
plot( x.axis, m5a$coef*100, ylim=c(-5, 15), pch=19, xlab="Bandwidth (days)",ylab="",main="Difference in Means", cex=.5, cex.axis=.9, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5a$lb*100, x1=x.axis, y1=m5a$ub*100)
axis(2, at=seq(-10,20,by=5), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)


plot( x.axis, m1a$coef*100, ylim=c(-5, 15), pch=19, xlab="Bandwidth (days)",ylab="",main="Linear", cex=.5, cex.axis=.8, cex.main=.9, yaxt="n", xaxt="n")
segments(x0=x.axis, y0=m1a$lb*100, x1=x.axis, y1=m1a$ub*100)
abline(h=0, col="red")
axis(2, at=seq(-10,20,by=5), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)

plot( x.axis, m2a$coef*100, ylim=c(-5, 15), pch=19, xlab="Bandwidth (days)",ylab="",main="Quadratic", cex=.5, cex.axis=.9, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2a$lb*100, x1=x.axis, y1=m2a$ub*100)
axis(2, at=seq(-10,20,by=5), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)


plot( x.axis, m6a$coef*100, ylim=c(-5, 15), pch=19, xlab="Bandwidth (days)",ylab="",main="Cubic", cex=.5, cex.axis=.8, cex.main=1, yaxt="n", xaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6a$lb*100, x1=x.axis, y1=m6a$ub*100)
axis(2, at=seq(-10,20,by=5), las=2, cex.axis=.8)
mtext(2, line=1.9, text="Effect\n(Percentage Points)", cex=.7)
axis(1, at=seq(-10,30,by=5), las=1, cex.axis=.7)



dev.off()








#############################################
#############################################
#############################################
#############################################
############################################# ONLINE APPENDIX
#############################################
#############################################
#############################################
#############################################


##Table  A1: discont in prob of finding a weapon, all stops, 2008-2015

#source("batch/full_models_all_data_batch.R")##compute HAC var covar matrix

res<-full.data.models(data=s, dv="found.weapon",  short.data=F, vc="output/all_data_HAC_vc.Rdata", gen.vc=F)
##print LaTex table
res[["table"]]

res

##local estimates of same quantity - Figure A4
stops.new.short<-s[s$dayno>=(memo-30) & s$dayno<=(memo+30), ]
dim(stops.new.short)
#out<-models(data=stops.new.short, N=1:30, dv="found.weapon",model="ols",memo=memo, int="no", cluster=F)
#save(out, file="output/local_test_findweapon_allstops.Rdata")

load("output/local_test_findweapon_allstops.Rdata")


m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]




x.axis<-1:30
file.name<-"output/bandwidth_findweapon_allstops.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)

plot( x.axis, m6$coef*100, ylim=c(-2, 2), pch=19, xlab="Bandwidth (days)",ylab="Effect (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
axis(2, seq(-5,5,1), las=2, cex.axis=.8)


dev.off()






###############
############### Figure A3
###############

s$cpw<-I(s$suspected.crime=="cpw")
range(s$dayno)
out.hit<-loess.time.series(data=s, dv="cpw",group="dayno", fun=mean, cut=memo)




##retrieve data for plot 
means.all.day<-out.hit$means.all.day
dd1<-out.hit$dd1
dd2<-out.hit$dd2
fit1<-out.hit$lower$fit1
fit2<-out.hit$upper$fit2
i1<-out.hit$lower$i1
i2<-out.hit$upper$i2

model1<-predict(loess(cpw.fun ~ dayno, data=means.all.day[means.all.day$dayno<memo,], span=.3, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1h<-model1$fit
lb1h<-model1$fit-1.96*model1$s
ub1h<-model1$fit+1.96*model1$s
i1h <- order(means.all.day$dayno)




means.all.day$cols<-NA
means.all.day$cols[as.numeric(means.all.day$dayno)<memo]<-"grey"
means.all.day$cols[as.numeric(means.all.day$dayno)>=memo]<-"grey"

file.name<-"output/hit.rate_cpw.pdf"
pdf(file.name)
plot(means.all.day$dayno, means.all.day[,2], pch=19, col=means.all.day$cols, cex=.3, xlim=c(min(means.all.day$dayno, na.rm=T), max(means.all.day$dayno, na.rm=T)), xlab="", ylab="", main="Share of Stops Where Weapon is Suspected\nOver time", axes=F, cex.main=1.3, ylim=c(0,.42))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=T, las=2, col.axis="black"
)
mtext(side=2, "Daily Share of Stops w/ CPW", line=3)
mtext(side=1, "Date", line=3, cex=.8)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
text( x=mdy.date(month=3, year=2013, day=5)-440, y=.39,labels="Day of\nProcedural\nChange" )

lines(dd1$dayno[i1], fit1[i1], col="red", lwd=2)
lines(dd2$dayno[i2], fit2[i2], col="red", lwd=2)
arrows(mdy.date(month=3, year=2013, day=5)-400, .35,mdy.date(month=3, year=2013, day=5), .35, length=.13 )
dev.off()





######################
######################
######################
#### Numerator/Denominator Tests Using All Data (Tables A2 and A3)
######################
######################
######################


sw2<-na.omit(sw[,c("period","distance","lag.hit","dayno","weapon.hit","weekday","year","month")])
sw2$count<-1
head(sw2)

##total weapon stops
stops<-summaryBy(count~dayno, data=sw2, FUN=sum, na.rm=T)
##total weapon stops producing a weapon
hits<-summaryBy(weapon.hit~dayno, data=sw2, FUN=sum, na.rm=T)

sw2$period2<-NA
sw2$period2[sw2$period=="pre-memo"]<-0
sw2$period2[sw2$period=="post-memo"]<-1
table(sw2$period2)

temp<-sw2[,c("period2","weekday","year","month","lag.hit","distance","dayno")]
temp<-unique(temp)
dim(temp)

sw3<-merge(hits, stops, by="dayno",all.x=T)
sw3<-merge(sw3, temp, by="dayno",all.x=T)
colnames(sw3)<-c("dayno","weapon.hit.sum","count.sum","period","weekday","year","month","lag.hit","distance")
head(sw3)
dim(sw3)

sw3$period[sw3$period==0]<-"pre-memo"
sw3$period[sw3$period==1]<-"post-memo"
sw3$period<-factor(sw3$period, levels=c("pre-memo","post-memo"))
table(sw3$period)

head(sw3)

sw3<-sw3[order(sw3$dayno),]

##all data
res.num<-full.data.models(data=sw3, dv="weapon.hit.sum", short.data=F, short.data.cut=NULL,  vc=NULL, gen.vc=T)
##print LaTex table
res.num[["table"]]

res.denom<-full.data.models(data=sw3, dv="count.sum", short.data=F, short.data.cut=NULL,  vc=NULL, gen.vc=T)
##print LaTex table
res.denom[["table"]]


######################################
######################################
######################################
###################################### RESULTS w/ OPTIMUM BANDWIDTH
######################################
######################################
######################################


### Estimates using rdrobust function with automatically selected bandwidth.


##use all data
#m<-rdrobust(x=sw$distance, y=sw$weapon.hit)


# band<-rdbwselect(x=sw$distance, y=sw$weapon.hit)
# print(band)

#save(m, file="output/rdrobust.results.Rdata")
load("output/rdrobust.results.Rdata")
names(m)
print(m)
summary(m)

#######################
#######################
####Clustered SEs - Table A4
#######################
#######################


##Table A4


dv<-"weapon.hit"
data1<-na.omit(sw[,c("period","distance","weapon.hit","precinct")])
data2<-na.omit(sw[,c("period","distance","weapon.hit","precinct","year","month","weekday","lag.hit")])


se1<-rep(NA, 8)
names.models<-c("diff.means","diff.means+controls","linear","linear+controls","squared","squared+controls","cubic","cubic+controls")
names(se1)<-names.models

##diff in means

m1<-lm(data1[,dv]~period, data=data1)
se1[1]<-summary(m1)$coefficients["periodpost-memo","Std. Error"]

m1a<-lm(data2[,dv]~period+factor(year)+factor(month)+factor(weekday)+lag.hit, data=data2)
se1[2]<-summary(m1a)$coefficients["periodpost-memo","Std. Error"]


##local linear
m2<-lm(data1[,dv]~period+distance+period*distance, data=data1)
se1[3]<-summary(m2)$coefficients["periodpost-memo","Std. Error"]
m2a<-lm(data2[,dv]~period+distance+period*distance+factor(year)+factor(month)+factor(weekday)+lag.hit, data=data2)
se1[4]<-summary(m2a)$coefficients["periodpost-memo","Std. Error"]

##second order poly
m3<-lm(data1[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2), data=data1)
se1[5]<-summary(m3)$coefficients["periodpost-memo","Std. Error"]
m3a<-lm(data2[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+factor(year)+factor(month)+factor(weekday)+lag.hit, data=data2)
se1[6]<-summary(m3a)$coefficients["periodpost-memo","Std. Error"]


##cubic
m4<-lm(data1[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3), data=data1)
se1[7]<-summary(m4)$coefficients["periodpost-memo","Std. Error"]

m4a<-lm(data2[,dv]~period+distance+period*distance+I(distance^2)+period*I(distance^2)+I(distance^3)+period*I(distance^3)+factor(year)+factor(month)+factor(weekday)+lag.hit, data=data2)
se1[8]<-summary(m4a)$coefficients["periodpost-memo","Std. Error"]


models2<-list(m1, m1a, m2, m2a, m3, m3a, m4, m4a)


coefs<-NA
se2<-NA
vc2<-list(NA)
for(i in 1:length(models2)){
	coefs[i]<-models2[[i]]$coefficients["periodpost-memo"]
	if(i %in% c(1,3,5,7)){
	vc2[[i]]<-vcovCluster(models2[[i]], data1$precinct)
	se2[i]<-sqrt(vc2[[i]]["periodpost-memo","periodpost-memo"])
	models2[[i]]$se<-se2[i]}
	if(i %in% c(2,4,6,8)){
	vc2[[i]]<-vcovCluster(models2[[i]], data2$precinct)
	se2[i]<-sqrt(vc2[[i]]["periodpost-memo","periodpost-memo"])
	models2[[i]]$se<-se2[i]}
}

for(i in 1:length(coefs)){
	lb<-coefs[[i]]-1.96*se2[[i]]
	print(lb)
}
##no estimate loses significance

apsrtable(models2[[1]],models2[[2]], models2[[3]], models2[[4]], models2[[5]] , models2[[6]], models2[[7]], models2[[8]],  model.names=c("Difference in Means", "Difference in Means","Linear","Linear","Squared","Squared","Cubic","Cubic"),omitcoef=c(1,3:length(m4a$coefficients) ),coef.names=c("Change in hit rate"), digits=4, se="robust"   )




##local estimates
stops.new.short<-sw[sw$dayno>=(memo-30) & sw$dayno<=(memo+30), ]
data<-stops.new.short
dim(stops.new.short)
out<-models(data=stops.new.short, N=1:30, dv="weapon.hit",model="ols",memo=memo, int="no", cluster=T, clust.var="precinct")

m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

##view results
out[["m1"]]

x.axis<-1:30
res<-list(m1, m2, m3, m4, m5 ,m6)

## Figure A5

##All models
file.name<-"output/bandwidth_precinct_cluster.pdf"
pdf( file.name, width=6, height=6)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1))
plot( x.axis, m5$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m1$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m2$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m3$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m4$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic 
+ Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)

plot( x.axis, m6$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
axis(2, at=seq(-100,100,by=2), cex.axis=.8, las=2)


dev.off()






######################
######################
######################
### Other DVs
######################
######################
######################

##Table B1 and Figure B1

##arrest rate
#global


#source("batch/arrested_models.R")

res<-full.data.models(data=s, dv="arrested", short.data=F ,  vc="output/arrested_HAC_vc.Rdata")
##print LaTex table
res[["table"]]

#local
#out<-models(data=s, N=1:30, dv="arrested", memo=memo, placebo=F, model="ols", int="no", cluster=F)
#save(out, file="output/arrested_bandwidth.Rdata")
load("output/arrested_bandwidth.Rdata")
##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]


x.axis<-seq(1, length(m1$coef), by=1)


file.name<-"output/bandwidth_arrest.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


dev.off()



##contraband rate 

##Table B2 and Figure B2


#source("batch/contraband_models.R")

res<-full.data.models(data=s, dv="found.contraband",  short.data=F,  vc="output/contraband_HAC_vc.Rdata")
##print LaTex table
res[["table"]]


#local
#out<-models(data=s, N=1:30, dv="found.contraband", memo=memo, placebo=F, model="ols", int="no", cluster=F)
#save(out, file="output/contraband_bandwidth.Rdata")
load("output/contraband_bandwidth.Rdata")
##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

file.name<-"output/bandwidth_contraband.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


dev.off()



##summons rate
##Table B3 and Figure B3

#source("batch/summons_models.R)

res<-full.data.models(data=s, dv="summons.issued",  short.data=F,  vc="output/summons_HAC_vc.Rdata")
##print LaTex table
res[["table"]]



#local
#out<-models(data=s, N=1:30, dv="summons.issued", memo=memo, placebo=F, model="ols", int="no", cluster=F)
#save(out, file="output/summons_bandwidth.Rdata")
load("output/summons_bandwidth.Rdata")
##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

file.name<-"output/bandwidth_summons.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-5, 10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


dev.off()







#################
#################
#### HOMICIDES - Table B4
#################
#################



d<-read.csv("data/nypd_weekly_homicides_2011_2013.csv")
head(d)
d$num<-1:nrow(d)
d
treat.location<-d$num[d$date=="3/4/2013-3/10/2013"]
d$treat<-0
d$treat[d$num>= treat.location]<-1
table(d$treat)
d$weeks.since<-(1-treat.location):((nrow(d)-treat.location))
##check coding of distance variable
d$weeks.since[d$num==treat.location-1]#is this -1? yes
d$weeks.since[d$num==treat.location]#is this zero? yes
d$weeks.since[d$num==treat.location+1]#is this 1? yes
##coding correct

d$date<-as.character(d$date)
temp<-strsplit(d$date, "-")
d$start.day<-NA
d$start.month<-NA
d$start.year<-NA

d$stop.day<-NA
d$stop.month<-NA
d$stop.year<-NA

for(i in 1:length(temp)){
	
d$start.month[i]<-strsplit(temp[[i]][1], "/")[[1]][1]
d$start.day[i]<-strsplit(temp[[i]][1], "/")[[1]][2]
d$start.year[i]<-strsplit(temp[[i]][1], "/")[[1]][3]

d$stop.month[i]<-strsplit(temp[[i]][2], "/")[[1]][1]
d$stop.day[i]<-strsplit(temp[[i]][2], "/")[[1]][2]
d$stop.year[i]<-strsplit(temp[[i]][2], "/")[[1]][3]
	
}
d$start.year[d$start.year=="3013"]<-"2013"##fix typo

head(d)



m1<-lm(homicides~treat, data=d)
m2<-lm(homicides ~treat+weeks.since+treat*weeks.since, data=d)
m3<-lm(homicides ~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2), data=d)
m4<-lm(homicides ~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+I(weeks.since^3)+treat*I(weeks.since^3), data=d)


m1a<-lm(homicides ~treat+factor(start.month), data=d)
m2a<-lm(homicides ~treat+weeks.since+treat*weeks.since+factor(start.month), data=d)
m3a<-lm(homicides ~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+factor(start.month), data=d)
m4a<-lm(homicides ~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+I(weeks.since^3)+treat*I(weeks.since^3)+factor(start.month), data=d)

apsrtable(m1,m1a,  m2, m2a, m3,m3a,  m4,m4a,  model.names=c("Difference in Means", "Difference in Means","Linear","Linear","Squared","Squared","Cubic","Cubic"),omitcoef=c(1,3:length(m4a$coefficients) ),coef.names=c("Discontinuity in Homicide Rate")   )



###ROBBERIES (for appendix)

##fit pre-treat loess for robberies
model1<-predict(loess(robberies ~ weeks.since, data=d[d$treat==0,], span=.5, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1h<-model1$fit
lb1h<-model1$fit-1.96*model1$s
ub1h<-model1$fit+1.96*model1$s
i1h <- order(d[d$treat==0,"weeks.since"])


##fit post-treat loess
model2<-predict(loess(robberies ~ weeks.since, data=d[d$treat==1,], span=.8, model=TRUE, x=TRUE, y=TRUE), se=T)
fit2h<-model2$fit
lb2h<-model2$fit-1.96*model2$s
ub2h<-model2$fit+1.96*model2$s
i2h <- order(d[d$treat==1,"weeks.since"])


#factor(start.month)*factor(start.year)
d$my<-(m<-lm(robberies~factor(start.month), data=d))$residuals
summary(lm(robberies~weeks.since+factor(start.month), data=d))

##fit pre-treat loess
model1<-predict(loess(my ~ weeks.since, data=d[d$treat==0,], span=.5, model=TRUE, x=TRUE, y=TRUE), se=T)
fit1hh<-model1$fit
lb1hh<-model1$fit-1.96*model1$s
ub1hh<-model1$fit+1.96*model1$s
i1hh <- order(d[d$treat==0,"weeks.since"])

##fit post-treat loess
model2<-predict(loess(my ~ weeks.since, data=d[d$treat==1,], span=.8, model=TRUE, x=TRUE, y=TRUE), se=T)
fit2hh<-model2$fit
lb2hh<-model2$fit-1.96*model2$s
ub2hh<-model2$fit+1.96*model2$s
i2hh <- order(d[d$treat==1,"weeks.since"])

##Figure B4


file.name<-"output/robbery.pdf"
pdf(file.name, width=12, height=5)

par(mar=c(4,6,2,4), mfrow=c(1,2))
plot(d$weeks.since, d$robberies, col="grey",cex=.6, axes=F,pch=19, ylab="Weekly Robberies",xlab="Weeks From Intervention", main="Weekly Robberies", cex.main=1, cex.axis=.5)
lines(d$weeks.since[d$treat==0][i1h], fit1h[i1h], col="black",lwd=2)
lines(d$weeks.since[d$treat==1][i2h], fit2h[i2h], col="black",lwd=2)
abline(v=0, lty=2, col="black", lwd=1)
axis(1, at=seq(-200, 100, by=50) , cex.axis=1)
axis(2, at=seq(0,500,100),las=2, cex.axis=1)
#abline(v=d$weeks.since[d$start.month%in%c(6,7,8)], col="grey",lty=2)
#axis(3, at=d$weeks.since,labels=d$month.lab, tck=0, cex.axis=.2)


plot(d$weeks.since, d$my, col="grey",cex=.6, axes=F,pch=19, ylab="Weekly Robberies\n(Controlling for Month)",xlab="Weeks From Intervention", main="Weekly Robberies\n(Controlling for Month)", cex.main=1, cex.axis=.5)
lines(d$weeks.since[d$treat==0][i1h], fit1hh[i1hh], col="black",lwd=2)
lines(d$weeks.since[d$treat==1][i2h], fit2hh[i2hh], col="black",lwd=2)
abline(v=0, lty=2, col="black", lwd=1)
axis(1, at=seq(-200, 100, by=50) , cex.axis=1)
axis(2, at=seq(-500,500,100),las=2, cex.axis=1)
dev.off()



m1<-lm(robberies~treat, data=d)
m2<-lm(robberies~treat+weeks.since+treat*weeks.since, data=d)
m3<-lm(robberies~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2), data=d)
m4<-lm(robberies~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+I(weeks.since^3)+treat*I(weeks.since^3), data=d)


m1a<-lm(robberies~treat+factor(start.month), data=d)
m2a<-lm(robberies~treat+weeks.since+treat*weeks.since+factor(start.month), data=d)
m3a<-lm(robberies~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+factor(start.month), data=d)
m4a<-lm(robberies~treat+weeks.since+treat*weeks.since+I(weeks.since^2)+treat*I(weeks.since^2)+I(weeks.since^3)+treat*I(weeks.since^3)+factor(start.month), data=d)


apsrtable(m1,m1a,  m2, m2a, m3,m3a,  m4,m4a,  model.names=c("Difference in Means", "Difference in Means","Linear","Linear","Squared","Squared","Cubic","Cubic"),omitcoef=c(1,3:length(m4a$coefficients) ),coef.names=c("Discontinuity in Robbery Rate")   )





###############################
###############################
###############################
### Change in Racial Composition of Stops + furtive movements
###############################
###############################
###############################




s$white2<-I(s$race2=="white")
s$black2<-I(s$race2=="black")
s$hisp2<-I(s$race2=="hispanic")
s$asian2<-I(s$race2=="asian")


table(sw$race2)

sw$white2<-I(sw$race2=="white")
sw$black2<-I(sw$race2=="black")
sw$hisp2<-I(sw$race2=="hispanic")
sw$asian2<-I(sw$race2=="asian")



##Figure B5

file.name<-"output/race_stops_all.pdf"
pdf(file.name)

par(mar=c(4,4.5,4,4), mfrow=c(3,2))


out.hit<-loess.time.series(data=s, dv="black2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,  xlab="", ylab="", main="Black Suspects", axes=F, cex.main=1  ,ylim=c(.4,.6))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=s, dv="white2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="White Suspects", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=s, dv="hisp2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Hispanic Suspects", axes=F, cex.main=1 , ylim=c(.2,.4))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)

out.hit<-loess.time.series(data=s, dv="asian2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Asian Suspects", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=s, dv="stopped.bc.furtive",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Furtive Movements", axes=F, cex.main=1 , ylim=c(.2,.6))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)

out.hit<-loess.time.series(data=s, dv="stopped.bc.object",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Suspicious Object", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


dev.off()




##Figure B6


file.name<-"output/race_stops_sw.pdf"
pdf(file.name)

par(mar=c(4,4.5,4,4), mfrow=c(3,2))


out.hit<-loess.time.series(data=sw, dv="black2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,  xlab="", ylab="", main="Black Suspects", axes=F, cex.main=1  ,ylim=c(.5,.8))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=sw, dv="white2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="White Suspects", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=sw, dv="hisp2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Hispanic Suspects", axes=F, cex.main=1 , ylim=c(.2,.4))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)

out.hit<-loess.time.series(data=sw, dv="asian2",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Asian Suspects", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


out.hit<-loess.time.series(data=sw, dv="stopped.bc.furtive",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Furtive Movements", axes=F, cex.main=1 , ylim=c(.35,.8))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)

out.hit<-loess.time.series(data=sw, dv="stopped.bc.object",group="dayno", fun=mean, cut=memo)
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="grey", cex=.3,   xlab="", ylab="", main="Suspicious Object", axes=F, cex.main=1 , ylim=c(0,.2))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.1), labels=paste(seq(0,100,by=10),"%",sep=""), las=2, col.axis="black")
mtext(side=2, "Percent of Stops", line=3, cex=.7)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="red", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="red", lwd=2)


dev.off()







###############################
###############################
###############################
###############################
############################### HETEROGENOUS TREATMENT EFFECTS - Tables C1-C3, Figures C1-C3
###############################
###############################
###############################
###############################

sb<-subset(sb, sb$suspected.crime=="cpw")
dim(sb)

sb<-sb[order(sb$dayno),]

sb$wpct<-as.numeric(as.character(sb$white))/as.numeric(as.character(sb$total))
summary(sb$wpct)

load("data/med.white.Rdata")##load median % white bg 
cut<-med.white


memo<-mdy.date(year=2013, month=3, day=5)

sb$b<-NA
sb$b[sb$wpct<cut]<-0
sb$b[sb$wpct>=cut]<-1 ##these are very white places
table(sb$b)

sb<-sb[order(sb$dayno),]



sb<-sb[order(sb$dayno),]



sb$weekday<-weekdays(as.Date(sb$date))


sb$period<-NA
sb$period[sb$dayno>=memo]<-"post-memo"
sb$period[sb$dayno<memo]<-"pre-memo"
sb$period<-factor(sb$period, levels=c("pre-memo","post-memo"))
table(sb$period)

sb$distance<-sb$dayno-memo

sw2<-sb

sw2$int<-sw2$b
sw2<-sw2[order(sw2$dayno),]



## High vs. Low % White in Block Group
#source("batch/full_white_bg_batch.R")##compute HAC var covar matrix

res<-full.data.models(data=NULL, dv=NULL, int=T, short.data=F,  model.data="output/std_models_white_bg.Rdata",  vc="output/white_bg_HAC_vc.Rdata", gen.vc=F)
res[["table"]]


##local effects
#out<-models(data=sw2, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols", int="yes",int.var="b", cluster=F)
#save(out, file="output/white_bg_interact_bandwidth.Rdata")

load("output/white_bg_interact_bandwidth.Rdata")

##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]



x.axis<-seq(1, length(m1$coef), by=1)


file.name<-"output/bandwidth_white_bg.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


dev.off()







## High vs. Low Homicides per Capita in Precinct


##### heterogeneity by precinct level crime
source("batch/nyc_precinct_crime.R") ##load precinct level mean crime rates, according to two measures
head(d)
table(d$pct)
length(unique(d$pct))

##load Rosenfeld and Fornago data to get precinct population estimates
pop<-read.dta("data/NYC Precincts Long 112712.dta")
head(pop)
pop<-pop[,c("pct","year","pop")]
pop<-pop[pop$year==2010,]##get 2010 population
head(pop)

table(pop$pct)
length(unique(pop$pct))

class(d$pct)
class(pop$pct)
pop$pct<-as.character(pop$pct)

pop.pct<-as.numeric(c(unique(pop$pct), NA))[order(as.numeric(c(unique(pop$pct), NA)))]
d.pct<-as.numeric(c(unique(d$pct)))[order(as.numeric(c(unique(d$pct))))]
sw.pct<-as.numeric(unique(sw$precinct))[order(as.numeric(unique(sw$precinct)))]
length(pop.pct)
length(d.pct)
length(sw.pct)

##note, the rosenfeld and fornago data exclude the central park precinct, probably because it has no "residents" so it cant be scaled by population. that's fine. footnote this.

d$pct[which(!(d$pct%in%pop$pct))]

head(pop)
pop<-pop[,c("pct","pop")]

d<-merge(d, pop, by="pct", all.x=T)

##make variables into low/high levels based on terciles
d$hom<-d$mean.homicide08.12/d$pop
cut.hom<-quantile(d$hom, probs=c(.333,.5,.666), na.rm=T)

d$hom3<-NA
d$hom3[d$hom<cut.hom[2]]<-0##above and below median
d$hom3[d$hom>=cut.hom[2] ]<-1
table(d$hom3)


head(sw)
class(sw$precinct)
sw$precinct<-as.character(sw$precinct)
table(sw$precinct)
dim(sw)
sw<-merge(sw, d, by.x="precinct", by.y="pct", all.x=T)
dim(sw)

##compute this with and without 77, 78 and 88, which were redrawn in late 2012 (update: tried this using the sw3 data frame and it didnt make much difference in the pattern of results)

sw2<-sw[,c("weapon.hit","period","mean.homicide08.12","dayno","distance","hom3","precinct","weekday","year","month","lag.hit")]


sw3<-na.omit(sw[!(sw$precinct%in%c("77","78","88")),c("weapon.hit","period","mean.homicide08.12","dayno","distance","hom3","precinct","weekday","year","month","lag.hit")])


##have estimated with and without these precincts. makes virtually no difference.

sw2<-sw2[order(sw2$dayno),]
identical(sw2$dayno, sw2$dayno[order(sw2$dayno)])

save(sw2,file="output/sw_hom.Rdata")##this includes the extra precincts




###TABLE of RESULTS FULL HOM DATA
#source("batch/full_batch_hom.R")##compute HAC var covar matrix

res<-full.data.models(data=NULL, dv=NULL, int=T,  short.data=F, model.data="output/std_models_hom.Rdata",  vc="output/hom_HAC_vc.Rdata", gen.vc=F)
res[["table"]]




##Local estimates



##interaction models
##local estimates of interaction term (diff in fiff conveying diff in treatment effect between low and high-homicide precincts)
#out<-models(data=sw2, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols", int="yes",int.var="hom3")
#save(out, file="output/homicide_interact_bandwidth.Rdata")

load("output/homicide_interact_bandwidth.Rdata")


##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]


x.axis<-seq(1, length(m6$coef), by=1)



file.name<-"output/bandwidth_hom.pdf"
pdf(file=file.name)

#par(mfrow=c(2,3))
layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2, cex.axis=.8)




dev.off()



## White vs. non-White Suspects
#source("batch/full_batch_white_sus.R")##compute HAC var covar matrix
res<-full.data.models(data=NULL, dv=NULL, int=T,  model.data="output/std_models_white_sus.Rdata",  vc="output/white_sus_HAC_vc.Rdata", gen.vc=F)
res[["table"]]


###local


sw2<-subset(sw, sw$dayno>=(memo-30) & sw$dayno<(memo+30))
dim(sw2)

#out<-models(data=sw2, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols", int="yes",int.var="white", cluster=F)
#save(out, file="output/white_interact_bandwidth.Rdata")

load("output/white_interact_bandwidth.Rdata")


##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-seq(1, length(m1$coef), by=1)
x.axis2<-seq(1, length(m1$coef), by=1)

file.name<-"output/bandwidth_white_diff.pdf"
pdf(file=file.name)


layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)


plot( x.axis2, m6$coef*100, ylim=c(-30, 50), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis2, y0=m6$lb*100, x1=x.axis2, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=10), las=2, cex.axis=.8)



dev.off()






#########################################
#########################################
#########################################
######################################### Reporting Bias Checks - Figures D1-D5
#########################################
#########################################
#########################################


##Non-weapon stops
stops.new<-s

stops.new$cpw<-I(stops.new$suspected.crime=="cpw")
 
 ##global

stops.new.short<-stops.new[stops.new$dayno>=(memo-30) & stops.new$dayno<(memo+30) & stops.new$cpw==F, ]

#out<-models(data=stops.new.short, N=1:30, dv="found.weapon",memo=memo, model="ols",int="no", cluster=F)
#save(out, file="output/figure3_data_noncpw.Rdata")

load("output/figure3_data_noncpw.Rdata")



##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-seq(1, length(m1$coef), by=1)

file.name<-"output/bandwidth_noncpw.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-5,5), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

dev.off()






### UNIFORM ANALYSIS
stops.new.short<-sw[sw$dayno>=(memo-30) & sw$dayno<(memo+30) & sw$officer.uniform==T, ]
#out<-models(data=stops.new.short, N=1:30, dv="weapon.hit", memo=memo, placebo=F, model="ols", cluster=F, int=F)
#save(out, file="output/uniform_data.Rdata")
load("output/uniform_data.Rdata")


##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-1:30




file.name<-"output/bandwidth_uniform.pdf"
pdf(file=file.name)

layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m1$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m2$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m3$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

plot( x.axis, m4$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)


plot( x.axis, m6$coef*100, ylim=c(-5,10), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100, 100, by=2), las=2, cex.axis=.8)

dev.off()




########THE RECLASSIFICATION HYPOTHESIS (FIGURE 14)

########PROPENSITY SCORE ANALYSIS
####USE PRE TREAT PARAMETERS TO PREDICT POST-TREATMENT PROB OF BEING LABELED A WEAPON STOP (among non weapon stops)

##calculate the propensity of being a cpw stop in the 10 days prior to the memo


### change in number of crimes
stops.new<-s
##make indicator for weapon stop
stops.new$cpw<-0
stops.new$cpw[stops.new$suspected.crime=="cpw"]<-1
stops.new$cpw[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$cpw)

stops.new$notcpw<-0
stops.new$notcpw[stops.new$suspected.crime!="cpw"]<-1
stops.new$notcpw[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$notcpw)

stops.new$robbery<-0
stops.new$robbery[stops.new$suspected.crime=="robbery"]<-1
stops.new$robbery[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$robbery)

stops.new$burglary<-0
stops.new$burglary[stops.new$suspected.crime=="burglary"]<-1
stops.new$burglary[is.na(stops.new$suspected.crime)]<-NA
table(stops.new$burglary)

vars<-c("suspect.height","suspect.weight","suspect.age","observation.period")

for(i in 1:length(vars)){
	
	stops.new[,paste(vars[i],"2",sep="")]<-as.numeric(as.character(stops.new[,vars[i]]))
	new<-paste(vars[i],"2",sep="")
	stops.new[,new]<-I(stops.new[, new]>=median(stops.new[, new], na.rm=T))	
}


summary(stops.new$suspect.height2)
summary(stops.new$suspect.weight2)
summary(stops.new$suspect.age2)
summary(stops.new$observation.period2)


file.name<-"output/prop_score2_predict.pdf"
pdf(file=file.name,width=6, height=6)


par(mfrow=c(2,2))

##one week bandwidth
##7-day bandwidth

m<-7

stops.new $post<-NA
stops.new $post[stops.new $dayno>=(memo-m) & stops.new $dayno<memo]<-0
stops.new $post[stops.new $dayno>=memo & stops.new $dayno<(memo+m) ]<-1
table(stops.new $post)



d<-na.omit(stops.new[,c("cpw","inside.outside","precinct", "location.housing", "observation.period2", "stopped.bc.object", "stopped.bc.desc", "stopped.bc.casing", "stopped.bc.lookout", "stopped.bc.clothing", "stopped.bc.drugs", "stopped.bc.furtive", "stopped.bc.violent", "stopped.bc.bulge", "stopped.bc.other",  "additional.proximity",  "additional.associating", "additional.direction", "additional.highcrime", "additional.time", "additional.sights", "additional.other", "suspect.sex", "suspect.race", "suspect.age2", "suspect.height2", "suspect.weight2", "suspect.hair","suspect.build","officer.uniform","weekday","hour","dayno","post")])
dim(d)
d$id<-1:(length(d[,1]))
j<-order(d$id[d$post==0])
k<-order(d$id[d$post==1])


psout <- glm(cpw~ inside.outside + factor(precinct) + factor(location.housing) +as.numeric(observation.period2) + stopped.bc.object   +stopped.bc.desc+         stopped.bc.casing+         stopped.bc.lookout+         stopped.bc.clothing+        stopped.bc.drugs+        
stopped.bc.furtive+         stopped.bc.violent+         stopped.bc.bulge+ stopped.bc.other  +additional.proximity          +additional.associating+ additional.direction+ additional.highcrime     
+ additional.time            +additional.sights +          additional.other   +suspect.sex+    suspect.race              
 +suspect.age2+                      +suspect.height2+            
 +suspect.weight2            +suspect.hair              +                suspect.build +officer.uniform+factor(weekday)+factor(hour),family=binomial,data=d[d$post==0, ]) 
summary(psout)
d$pre.prop<-NA
d$pre.prop[d$post==0 ][j]<-psout$fitted.values[j]
summary(d$pre.prop)


##fit a model using the post-treatment observations
prob2<-predict(psout, newdata=d[d$post==1, ], type="response")

d$post.prop<-NA
d$post.prop[d$post==1 ][k]<-prob2[k]

plot(density(d$pre.prop[d$cpw==0 & d$post==0]), lty=1, ylim=c(0,8),main="7-Day Bandwidth", xlab="Propensity to be Labeled a Weapon Stop", yaxt="n")
lines(density(d$post.prop[d$cpw==0 & d$post==1]), lty=2, col="red")
legend("topright",legend=c("Pre-Memo","Post-Memo"), lty=c(1,2),col=c("black","red"))
abline(v=mean(d$pre.prop[d$cpw==0 & d$post==0], na.rm=T))
abline(v=mean(d$post.prop[d$cpw==0 & d$post==1], na.rm=T), lty=2, col="red")
axis(2, seq(-10,10,by=2), las=2)


##two week bandwidth
m<-14

stops.new $post<-NA
stops.new $post[stops.new $dayno>=(memo-m) & stops.new $dayno<memo]<-0
stops.new $post[stops.new $dayno>=memo & stops.new $dayno<(memo+m) ]<-1
table(stops.new $post)



d<-na.omit(stops.new[,c("cpw","inside.outside","precinct", "location.housing", "observation.period2", "stopped.bc.object", "stopped.bc.desc", "stopped.bc.casing", "stopped.bc.lookout", "stopped.bc.clothing", "stopped.bc.drugs", "stopped.bc.furtive", "stopped.bc.violent", "stopped.bc.bulge", "stopped.bc.other",  "additional.proximity",  "additional.associating", "additional.direction", "additional.highcrime", "additional.time", "additional.sights", "additional.other", "suspect.sex", "suspect.race", "suspect.age2", "suspect.height2", "suspect.weight2", "suspect.hair","suspect.build","officer.uniform","weekday","hour","dayno","post")])
dim(d)
d$id<-1:(length(d[,1]))
j<-order(d$id[d$post==0])
k<-order(d$id[d$post==1])


psout <- glm(cpw~ inside.outside + factor(precinct) + factor(location.housing) +as.numeric(observation.period2) + stopped.bc.object   +stopped.bc.desc+         stopped.bc.casing+         stopped.bc.lookout+         stopped.bc.clothing+        stopped.bc.drugs+        
stopped.bc.furtive+         stopped.bc.violent+         stopped.bc.bulge+ stopped.bc.other  +additional.proximity          +additional.associating+ additional.direction+ additional.highcrime     
+ additional.time            +additional.sights +          additional.other   +suspect.sex+    suspect.race              
 +suspect.age2+                      +suspect.height2+            
 +suspect.weight2            +suspect.hair              +                suspect.build +officer.uniform+factor(weekday)+factor(hour),family=binomial,data=d[d$post==0, ]) 
summary(psout)
d$pre.prop<-NA
d$pre.prop[d$post==0 ][j]<-psout$fitted.values[j]
summary(d$pre.prop)


##fit a model using the post-treatment observations
prob2<-predict(psout, newdata=d[d$post==1, ], type="response")

d$post.prop<-NA
d$post.prop[d$post==1 ][k]<-prob2[k]

plot(density(d$pre.prop[d$cpw==0 & d$post==0]), lty=1, ylim=c(0,8),main="14-Day Bandwidth", xlab="Propensity to be Labeled a Weapon Stop", yaxt="n")
lines(density(d$post.prop[d$cpw==0 & d$post==1]), lty=2, col="red")
legend("topright",legend=c("Pre-Memo","Post-Memo"), lty=c(1,2),col=c("black","red"))
abline(v=mean(d$pre.prop[d$cpw==0 & d$post==0], na.rm=T))
abline(v=mean(d$post.prop[d$cpw==0 & d$post==1], na.rm=T), lty=2, col="red")
axis(2, seq(-10,10,by=2), las=2)




##three week bandwidth

m<-21

stops.new $post<-NA
stops.new $post[stops.new $dayno>=(memo-m) & stops.new $dayno<memo]<-0
stops.new $post[stops.new $dayno>=memo & stops.new $dayno<(memo+m) ]<-1
table(stops.new $post)



d<-na.omit(stops.new[,c("cpw","inside.outside","precinct", "location.housing", "observation.period2", "stopped.bc.object", "stopped.bc.desc", "stopped.bc.casing", "stopped.bc.lookout", "stopped.bc.clothing", "stopped.bc.drugs", "stopped.bc.furtive", "stopped.bc.violent", "stopped.bc.bulge", "stopped.bc.other",  "additional.proximity",  "additional.associating", "additional.direction", "additional.highcrime", "additional.time", "additional.sights", "additional.other", "suspect.sex", "suspect.race", "suspect.age2", "suspect.height2", "suspect.weight2", "suspect.hair","suspect.build","officer.uniform","weekday","hour","dayno","post")])
dim(d)
d$id<-1:(length(d[,1]))
j<-order(d$id[d$post==0])
k<-order(d$id[d$post==1])


psout <- glm(cpw~ inside.outside + factor(precinct) + factor(location.housing) +as.numeric(observation.period2) + stopped.bc.object   +stopped.bc.desc+         stopped.bc.casing+         stopped.bc.lookout+         stopped.bc.clothing+        stopped.bc.drugs+        
stopped.bc.furtive+         stopped.bc.violent+         stopped.bc.bulge+ stopped.bc.other  +additional.proximity          +additional.associating+ additional.direction+ additional.highcrime     
+ additional.time            +additional.sights +          additional.other   +suspect.sex+    suspect.race              
 +suspect.age2+                      +suspect.height2+            
 +suspect.weight2            +suspect.hair              +                suspect.build +officer.uniform+factor(weekday)+factor(hour),family=binomial,data=d[d$post==0, ]) 
summary(psout)
d$pre.prop<-NA
d$pre.prop[d$post==0 ][j]<-psout$fitted.values[j]
summary(d$pre.prop)


##fit a model using the post-treatment observations
prob2<-predict(psout, newdata=d[d$post==1, ], type="response")

d$post.prop<-NA
d$post.prop[d$post==1 ][k]<-prob2[k]

plot(density(d$pre.prop[d$cpw==0 & d$post==0]), lty=1, ylim=c(0,8),main="21-Day Bandwidth", xlab="Propensity to be Labeled a Weapon Stop", yaxt="n")
lines(density(d$post.prop[d$cpw==0 & d$post==1]), lty=2, col="red")
legend("topright",legend=c("Pre-Memo","Post-Memo"), lty=c(1,2),col=c("black","red"))
abline(v=mean(d$pre.prop[d$cpw==0 & d$post==0], na.rm=T))
abline(v=mean(d$post.prop[d$cpw==0 & d$post==1], na.rm=T), lty=2, col="red")
axis(2, seq(-10,10,by=2), las=2)



##four week bandwidth
m<-30

stops.new $post<-NA
stops.new $post[stops.new $dayno>=(memo-m) & stops.new $dayno<memo]<-0
stops.new $post[stops.new $dayno>=memo & stops.new $dayno<(memo+m) ]<-1
table(stops.new $post)



d<-na.omit(stops.new[,c("cpw","inside.outside","precinct", "location.housing", "observation.period2", "stopped.bc.object", "stopped.bc.desc", "stopped.bc.casing", "stopped.bc.lookout", "stopped.bc.clothing", "stopped.bc.drugs", "stopped.bc.furtive", "stopped.bc.violent", "stopped.bc.bulge", "stopped.bc.other",  "additional.proximity",  "additional.associating", "additional.direction", "additional.highcrime", "additional.time", "additional.sights", "additional.other", "suspect.sex", "suspect.race", "suspect.age2", "suspect.height2", "suspect.weight2", "suspect.hair","suspect.build","officer.uniform","weekday","hour","dayno","post")])
dim(d)
d$id<-1:(length(d[,1]))
j<-order(d$id[d$post==0])
k<-order(d$id[d$post==1])


psout <- glm(cpw~ inside.outside + factor(precinct) + factor(location.housing) +as.numeric(observation.period2) + stopped.bc.object   +stopped.bc.desc+         stopped.bc.casing+         stopped.bc.lookout+         stopped.bc.clothing+        stopped.bc.drugs+        
stopped.bc.furtive+         stopped.bc.violent+         stopped.bc.bulge+ stopped.bc.other  +additional.proximity          +additional.associating+ additional.direction+ additional.highcrime     
+ additional.time            +additional.sights +          additional.other   +suspect.sex+    suspect.race              
 +suspect.age2+                      +suspect.height2+            
 +suspect.weight2            +suspect.hair              +                suspect.build +officer.uniform+factor(weekday)+factor(hour),family=binomial,data=d[d$post==0, ]) 
summary(psout)
d$pre.prop<-NA
d$pre.prop[d$post==0 ][j]<-psout$fitted.values[j]
summary(d$pre.prop)


##fit a model using the post-treatment observations
prob2<-predict(psout, newdata=d[d$post==1, ], type="response")

d$post.prop<-NA
d$post.prop[d$post==1 ][k]<-prob2[k]

plot(density(d$pre.prop[d$cpw==0 & d$post==0]), lty=1, ylim=c(0,8),main="30-Day Bandwidth", xlab="Propensity to be Labeled a Weapon Stop", yaxt="n")
lines(density(d$post.prop[d$cpw==0 & d$post==1]), lty=2, col="red")
legend("topright",legend=c("Pre-Memo","Post-Memo"), lty=c(1,2),col=c("black","red"))
abline(v=mean(d$pre.prop[d$cpw==0 & d$post==0], na.rm=T))
abline(v=mean(d$post.prop[d$cpw==0 & d$post==1], na.rm=T), lty=2, col="red")
axis(2, seq(-10,10,by=2), las=2)


dev.off()


##tests referenced in Appendix D

##KS-Test using 30-day bandwidth

##Compute quantiles referenced in Appendix D

ks.test(d$pre.prop, d$post.prop)
quantile(d$pre.prop, probs=c(.25,.5,.75), na.rm=T)
quantile(d$post.prop, probs=c(.25,.5,.75), na.rm=T)







######################
######################
######################
### Rate of ID refusals
######################
######################
######################

s$refused<-I(s$identification=="refused")


out.hit<-loess.time.series(data=s, dv="refused",group="dayno", fun=mean, cut=memo)

##retrieve data for plot 1
means.all.day<-out.hit$means.all.day
dd1<-out.hit$dd1
dd2<-out.hit$dd2
fit1<-out.hit$lower$fit1
fit2<-out.hit$upper$fit2
i1<-out.hit$lower$i1
i2<-out.hit$upper$i2


upper<-quantile(out.hit$means.all.day$refused.fun, probs=c(.99))


#generate figure
file.name<-paste(wd, "output/refused_ID.pdf", sep="")
pdf(file.name)

plot(means.all.day$dayno, means.all.day[,2], pch=19, col="grey", cex=.3, xlim=c(min(means.all.day$dayno, na.rm=T), max(means.all.day$dayno, na.rm=T)), xlab="", ylab="", main="Rate of ID Refusals Over Time", axes=F, cex.main=1.3, ylim=c(0,upper))
axis(1, col="black", at=seq(17532, 20453, by=365*2), labels=dates,  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.01), labels=T, las=2, col.axis="black"
)
mtext(side=2, "Daily Rate of ID Refusals", line=3)
mtext(side=1, "Date", line=3, cex=.8)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
text( x=mdy.date(month=3, year=2013, day=5)-440, y=.49,labels="Day of\nProcedural\nChange" )
lines(dd1$dayno[i1], fit1[i1], col="black", lwd=2)
lines(dd2$dayno[i2], fit2[i2], col="black", lwd=2)
arrows(mdy.date(month=3, year=2013, day=5)-400, .43,mdy.date(month=3, year=2013, day=5), .43, length=.13 )

dev.off()




########################
########################
########################
########################
###Officers in uniform, all data, Figure D5
########################
########################
########################
########################

#source("batch/full_uniform.R")

load("output/std_models_uniform.Rdata")
load("output/uniform_HAC_vc.Rdata")

##Input point estimates from model using all weapon stops (Table 1)
coefs<-c(.051 , .031, .030, .022, .029, .020, .013, .010)*100
ses<-c(.002, .003, .002, .003, .003, .003, .003, .004)*100
lbs<-coefs-1.96*ses
ubs<-coefs+1.96*ses
model.names<-c("Diff. in Means","Diff. in Means\n + Controls","Linear","Linear\n + Controls", "Quadratic","Quadratic\n + Controls","Cubic","Cubic \n + Controls")

coefs2<-NA
ses2<-NA
lbs2<-NA
ubs2<-NA
for(i in 1:length(models2)){
coefs2[i]<-models2[[i]]$coefficients["periodpost-memo"]*100

temp<-summary(models2[[i]])$coefficients["periodpost-memo","Std. Error"]
temp2<-sqrt(vc[[i]]["periodpost-memo","periodpost-memo"])

if(temp>temp2){se<-temp}
if(temp<temp2){se<-temp2}


ses2[i]<-se*100
lbs2[i]<-coefs2[i]-1.96*ses2[i]
ubs2[i]<-coefs2[i]+1.96*ses2[i]

}

coefs.all<-c(coefs[1],coefs2[1],coefs[2],coefs2[2], coefs[3],coefs2[3],coefs[4],coefs2[4],coefs[5],coefs2[5],coefs[6],coefs2[6],coefs[7],coefs2[7],coefs[8],coefs2[8])


ses.all<-c(ses[1],ses2[1],ses[2],ses2[2], ses[3],ses2[3],ses[4],ses2[4],ses[5],ses2[5],ses[6],ses2[6],ses[7],ses2[7],ses[8],ses2[8])

lbs.all<-c(lbs[1],lbs2[1],lbs[2],lbs2[2], lbs[3],lbs2[3],lbs[4],lbs2[4],lbs[5],lbs2[5],lbs[6],lbs2[6],lbs[7],lbs2[7],lbs[8],lbs2[8])

ubs.all<-c(ubs[1],ubs2[1],ubs[2],ubs2[2], ubs[3],ubs2[3],ubs[4],ubs2[4],ubs[5],ubs2[5],ubs[6],ubs2[6],ubs[7],ubs2[7],ubs[8],ubs2[8])



model.names2<-c("Diff. in Means","","Diff. in Means\n + Controls","","Linear","","Linear\n + Controls", "","Quadratic","","Quadratic\n + Controls","","Cubic","","Cubic \n + Controls","")




file.name<-"output/all_models_coefs_uniform.pdf"
pdf(file.name, height=6, width=6)

shape<-c(19, 15)
par(mar=c(4,7,4,4))
y.axis<-length(coefs.all):1
plot(coefs.all, y.axis, pch=shape, xlim=c(0, max(ubs.all)), axes=F, xlab="Change in Hit Rate (Percentage Points)",ylab="",main="Change in Hit Rate on Day of Intervention\n (All Weapon Stops, 2008-2015)", col=c("black","blue"))
axis(1, at=seq(0, .1, .01)*100)
axis(1, at=seq(0, .1, .01)*100, labels=F, tck=0)
axis(1, at=seq(0, .1, .01)*100, labels=F, tck=0)

axis(2, at=y.axis-.2, labels=model.names2, las=2, tck=0)
axis(2, at=seq(-10,20,by=1), labels=F, las=2, tck=0)
axis(2, at=seq(-10,20,by=1), labels=F, las=2, tck=0)
axis(2, at=seq(-10,20,by=1), labels=F, las=2, tck=0)

abline(v=0, lty=2)
segments(lbs.all, y.axis, ubs.all, y.axis,col=c("black","blue"))
abline(h=seq(2.5, 15.5, by=2), col="grey")
legend("bottomright",legend=c("All officers","Officers in uniform"),pch=c(19,15),lty=1,bg="white",col=c("black","blue"))

dev.off()









##########################
##########################
##########################
########################## PLACEBO CHECKS - Appendix E
##########################
##########################
##########################




#######HYPOTHETICAL CUTOFFS in PRETREATMENT DATA


memo<-mdy.date(month=3, year=2013, day=5)#memo



#subset to every day prior to actual intervention that might be used as a cutoff 
stops.new.cpw<-sw[sw$suspected.crime=="cpw" & sw$dayno<=memo, ]
dim(stops.new.cpw)


#make storage vectors for coefficients
beta.loc.lin<-NA
beta.loc.lin.sq<-NA
beta.loc.lin.lag<-NA
beta.loc.lin.lag.sq<-NA
beta.cubic<-NA
beta.diff<-NA

results.day<-list(beta.loc.lin, beta.loc.lin.sq,beta.loc.lin.lag,beta.loc.lin.lag.sq, beta.cubic, beta.diff  )
names(results.day)<-c("beta.loc.linear", "beta.loc.linear.sq","beta.loc.linear.lag","beta.loc.linear.lag.sq", "beta.cubic", "beta.diff" )


day.lengths<-c(15, 30)
length(day.lengths)


results<-list(NA)
results[[1]]<-results.day


for(i in 1:length(day.lengths)){
	
	results[[i]]<-results.day
	
	
}

names(results)<-as.character(day.lengths)


##ESTIMATE PLACEBO EFFECTS - Figures E1 and E2
#source("batch/hyp_cutoff_batch.R")

##load observed effects using actual cutoff date
 
load("output/figure3_data.Rdata")

##view results
##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]


titles<-c("Difference\nin Means","Linear","Squared","Cubic")

file.name<-"output/placebo_dist_all_models.pdf"
pdf(file.name)

##LOAD PLACEBO EFFECTS
load("output/placebo_meta_new2.Rdata")

day.lengths<-c(15,30)
#compute observed effect



par(mfrow=c(2,4))
for(j in 1:length(day.lengths)){

##diff means

observed<-NA
observed[1]<-m5$coef[day.lengths[j]]##15 day bandwidth

 hist(results[[j]]$beta.diff*100, main=titles[1], xlab="Effect (Percentage Points)", cex.axis=.5, cex.lab=1, cex.main=1, xlim=c(-5,6), yaxt="n")
 abline(v=quantile(results[[j]]$beta.diff*100, probs=c(.025,.975)), lty=2)
 abline(v=observed[1]*100, col="red")
axis(2, at=seq(0, 500, by=50), las=2, cex.axis=.5)


# ##local linear
observed[2]<-m1$coef[day.lengths[j]]##15 day bandwidth
 hist(results[[j]]$beta.loc.linear*100, main=titles[2], xlab="Effect (Percentage Points)", cex.axis=.5, cex.lab=1, cex.main=1,, xlim=c(-5,6), yaxt="n")
 abline(v=quantile(results[[j]]$beta.loc.linear*100, probs=c(.025,.975)), lty=2)
 abline(v=observed[1]*100, col="red")
 axis(2, at=seq(0, 500, by=50), las=2, cex.axis=.5)


# ##local linear sq
observed[3]<-m2$coef[day.lengths[j]]##15 day bandwidth
 hist(results[[j]]$beta.loc.linear.sq*100, main=titles[3], xlab="Effect (Percentage Points)", cex.axis=.5, cex.lab=1, cex.main=1, xlim=c(-5,6), yaxt="n")
 abline(v=quantile(results[[j]]$beta.loc.linear.sq*100, probs=c(.025,.975)), lty=2)
 abline(v=observed[1]*100, col="red")
 axis(2, at=seq(0, 500, by=50), las=2, cex.axis=.5)

# ##cubic
observed[4]<-m6$coef[day.lengths[j]]##15 day bandwidth
 hist(results[[j]]$beta.cubic*100, main=titles[4], xlab="Effect (Percentage Points)", cex.axis=.5, cex.lab=1, cex.main=1, xlim=c(-5,6), yaxt="n")
 abline(v=quantile(results[[j]]$beta.cubic*100, probs=c(.025,.975)), lty=2)
 abline(v=observed[1]*100, col="red")
 axis(2, at=seq(0, 500, by=50), las=2, cex.axis=.5)



}

dev.off()






##REPEAT USING ALL DATA


# 
model.names<-c("diff","diff+controls","loc. linear","loc. linear+controls","squared", "squared+controls","cubic","cubic+controls")
results<-list(NA)

coefs<-NA
ses<-NA


days<-as.date(unique(sw$dayno))
days<-days[order(days)]
days<-days[days<as.date("4Mar2013")]##make march 3 the last cutoff, leaving march 4 as the "treatment" group
days<-days[days>as.date("2Jan2008")]##start at Jan 2 so lag.hit is not missing
head(days)
length(days)
summary(days)

for(i in 1:length(model.names)){
	results[[i]]<-cbind.data.frame(model.name=rep(model.names[i],length(days)), coef=rep(NA, length(days)), se=rep(NA, length(days)))
}


names(results)<-model.names


stops.new.short<-subset(sw, sw$dayno<(memo))
dim(stops.new.short)

#ESTIMATE PLACEBO EFFECTS
#source("batch/hyp_cutoff_batch2.R")

##LOAD PLACEBO EFFECTS
load("output/placebo_all.Rdata")

x.axis<-seq(0, length(m1$coef)-1, by=1)

##get observed discont using all data


observed<-NA

m<-lm(weapon.hit~period+distance+period*distance, data=sw)
observed[1]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[1]<-"loc. linear"

m<-lm(weapon.hit~period+distance+period*distance+factor(year)+factor(month)+factor(weekday)+lag.hit, data=sw)
observed[2]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[2]<-"loc. linear+controls"

m<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2), data=sw)
observed[3]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[3]<-"squared"


m<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2)+factor(year)+factor(month)+factor(weekday)+lag.hit, data=sw)
observed[4]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[4]<-"squared+controls"



m<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3), data=sw)
observed[5]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[5]<-"cubic"


m<-lm(weapon.hit~period+distance+period*distance + I(distance^2)+period*I(distance^2) +  I(distance^3)+period*I(distance^3)+factor(year)+factor(month)+factor(weekday)+lag.hit, data=sw)
observed[6]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[6]<-"cubic+controls"

m<-lm(weapon.hit~period, data=sw)
observed[7]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[7]<-"diff"

m<-lm(weapon.hit~period+factor(year)+factor(month)+factor(weekday)+lag.hit, data=sw)
observed[8]<-summary(m)$coefficients["periodpost-memo","Estimate"]
names(observed)[8]<-"diff+controls"

names(observed)

k<-c("diff" ,"diff+controls" ,"loc. linear" ,"loc. linear+controls" ,"squared" ,"squared+controls","cubic" ,               "cubic+controls")
main.titles<-c("Difference in Means" ,"Difference in Means\n+Controls" ,"Linear" ,"Linear+Controls" ,"Squared" ,"Squared+Controls","Cubic" , "Cubic+Controls")      



results<-results[k]
observed<-observed[k]

file.name<-"output/placebo_alldata.pdf"
pdf(file=file.name)

par(mfrow=c(2,4))

for(i in 1:length(names(observed))){
	
 hist(results[[names(observed)[i]]]$coef*100, main=main.titles[i], xlab="Effect (Percentage Points)", cex.axis=.5, cex.lab=1, cex.main=1, xlim=c(-5,6), yaxt="n", breaks=20)
 abline(v=observed[names(observed)[i]]*100, lty=1, col="red")
 abline(v=quantile(results[[names(observed)[i]]]$coef, probs=c(.025, .975), na.rm=T)*100, lty=2)
 axis(2, at=seq(0, 2000, by=200), las=2, cex.axis=.5)
}

	


dev.off()



########################
############
##### Anticipatory Behavior - Figure E3
############
############ Zoom in on Hit Rate Right Before Intervention

out.hit<-loess.time.series(data=sw[sw$dayno>=(memo-200)&sw$dayno<(memo)+200,], dv="weapon.hit",group="dayno", fun=mean, cut=memo)

##Figure E3

##zoom in on area around intervention
file.name<-paste(wd, "output/hit.rate.time_loess_zoom.pdf", sep="")
pdf(file=file.name)
par(mar=c(4,4.5,4,4))
plot(out.hit$means.all.day$dayno, out.hit$means.all.day[,2], pch=19, col="dimgrey", cex=.3,  xlab="", ylab="", main="Hit Rate Over Time", axes=F, cex.main=1.5, ylim=c(0,.15), xlim=c(as.numeric(memo)-200, as.numeric(memo)+200))
axis(1, col="black", at=seq(as.numeric(memo)-500, as.numeric(memo)+500, 50), labels=seq(-500, 500, 50),  cex.axis=.8)
axis(2, col="black", at=seq(0,1, by=.05), labels=T, las=2, col.axis="black")
mtext(side=2, "Daily Weapon Recovery Rate", line=3, cex=1.5)
mtext(side=1, "Days From Intervention", line=3, cex=1.5)
abline(v=mdy.date(month=3, year=2013, day=5), lwd=1,lty=2,  col="black")##march 5, 2013 - memo
text( x=mdy.date(month=3, year=2013, day=5)-440, y=.29,labels="Day of\nProcedural\nChange" )
lines(out.hit$dd1$dayno[out.hit$lower$i1], out.hit$lower$fit1[out.hit$lower$i1], col="black", lwd=2)
lines(out.hit$dd2$dayno[out.hit$upper$i2], out.hit$upper$fit2[out.hit$upper$i2], col="black", lwd=2)
arrows(mdy.date(month=3, year=2013, day=5)-400, .24,mdy.date(month=3, year=2013, day=5), .24 )
dev.off()




############# Placebo discontinuity using 60 days prior 

cut<-mdy.date(month=3, year=2013, day=5)#define intervention day, March 5, 2013

stops.new.short<-sw[sw$dayno>=(cut-60) & sw$dayno<(cut), ]

stops.new.short$period2<-I(stops.new.short$dayno>=(cut-30))
stops.new.short$period3<-NA
stops.new.short$period3[stops.new.short$dayno<(cut-30)]<-0
stops.new.short$period3[stops.new.short$dayno>=(cut-30)]<-1
table(stops.new.short$period3)
table(stops.new.short$period2)
cor(stops.new.short$period2, stops.new.short$period3)


length(table(stops.new.short$dayno))##60 days
length(unique(stops.new.short$dayno[stops.new.short$period2==T]))##30 days post
length(unique(stops.new.short$dayno[stops.new.short$period2==F]))##30 days pre
max(stops.new.short$dayno)##final day is day before treatment


##function revised to include both HAC and standard SEs
#out<-models(data=stops.new.short, N=1:30, dv="weapon.hit", model="ols",memo=(cut-30), placebo=T, int="no", cluster=F)###
#save(out, file="output/local_discont_placebo.Rdata")
load(file="output/local_discont_placebo.Rdata")
##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-seq(1, length(m1$coef), by=1)
res<-list(m1, m2, m3, m4, m5 ,m6)



##Figure E4

##All models
file.name<-"output/bandwidth_placebo_local.pdf"
pdf( file.name, width=6, height=6)
layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m1$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m2$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m3$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m4$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial 
+ Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m6$coef*100, ylim=c(-10, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
axis(2, at=seq(-100,100, by=2), las=2)

dev.off()







################# BRONX PLACEBO TEST - Table E1, Figure E5


##clean halls placebo check

bronx.pre<-as.character(c(40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 52))##source: https://en.wikipedia.org/wiki/Organization_of_the_New_York_City_Police_Department , Note: wikipedia retrieved these from nypd in 2012 so these are correct for bronx precincts at that time

sw$bronx<-NA
sw$bronx[sw$precinct%in%bronx.pre]<-1
sw$bronx[!(sw$precinct%in%bronx.pre)]<-0
table(sw$bronx)

memo2<-mdy.date(month=1, day=8, year=2013)##date of court ruling


sw2<-subset(sw, sw$dayno>=(memo2-30) & sw$dayno<(memo2+30) & sw$bronx==1)
dim(sw2)



 #out<-models(data=sw2, N=1:30, dv="weapon.hit", memo=memo2, placebo=T, model="ols", int="no", cluster=F)
 #save(out, file="output/bandwidth_bronx_placebo.Rdata")

load("output/bandwidth_bronx_placebo.Rdata")

##reorganize results for plotting
m1<-out[["m1"]]
m2<-out[["m2"]]
m3<-out[["m3"]]
m4<-out[["m4"]]
m5<-out[["m5"]]
m6<-out[["m6"]]

x.axis<-seq(1, length(m1$coef), by=1)


file.name<-"output/bandwidth_bronx_placebo.pdf"
pdf(file=file.name)

#par(mfrow=c(2,3))
layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE), widths=c(1,1,1), heights=c(1,1,.75))

plot( x.axis, m5$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Difference in Means", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m5$lb*100, x1=x.axis, y1=m5$ub*100)
points( x.axis, m5$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m1$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear", cex=.7, cex.axis=.8, yaxt="n")
segments(x0=x.axis, y0=m1$lb*100, x1=x.axis, y1=m1$ub*100)
abline(h=0, col="red")
points( x.axis, m1$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m2$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Second Order Polynomial", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m2$lb*100, x1=x.axis, y1=m2$ub*100)
points( x.axis, m2$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m3$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Local Linear +
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m3$lb*100, x1=x.axis, y1=m3$ub*100)
points( x.axis, m3$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)

plot( x.axis, m4$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Quadratic
Prior Day's Hit Rate", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m4$lb*100, x1=x.axis, y1=m4$ub*100)
points( x.axis, m4$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)


plot( x.axis, m6$coef*100, ylim=c(-15, 15), pch=19, xlab="Bandwidth (days)",ylab="Effect Size (Percentage Points)",main="Cubic", cex=.7, cex.axis=.8, yaxt="n")
abline(h=0, col="red")
segments(x0=x.axis, y0=m6$lb*100, x1=x.axis, y1=m6$ub*100)
points( x.axis, m6$coef*100, cex=.7)
axis(2, at=seq(-100,100, by=2), las=2)



dev.off()





###Global Results

bronx.pre<-as.character(c(40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 52))

sw$bronx<-NA
sw$bronx[sw$precinct%in%bronx.pre]<-1
sw$bronx[!(sw$precinct%in%bronx.pre)]<-0
table(sw$bronx)

memo2<-mdy.date(month=1, day=8, year=2013)
memo<-mdy.date(month=3, day=5, year=2013)



sw2<-subset(sw, sw$dayno<memo & sw$bronx==1 )##subset to days prior to actual intervention for precincts in the bronx

sw2$period<-NA
sw2$period[sw2$dayno>=memo2]<-"post-memo"
sw2$period[sw2$dayno<memo2]<-"pre-memo"
sw2$period<-factor(sw2$period,levels=c("pre-memo","post-memo"))
table(sw2$period)


sw2$distance<-NA
sw2$distance<-sw2$dayno-memo2

res<-full.data.models(data=sw2, int=F, gen.vc=F,short.data=F, dv="weapon.hit",   vc="output/bronx_HAC_vc.Rdata")
##print LaTex table
res[["table"]]











